module functions_mod
USE parameters_mod
USE parameters_energy_mod
IMPLICIT NONE
CONTAINS




FUNCTION c_upper_fnc2(labor,effective_wage,capital,kprime,kdim,tr,r)
INTEGER, intent(in)   ::kdim
DOUBLE PRECISION,    intent(in)   ::kprime(kdim),capital,labor,effective_wage
DOUBLE PRECISION,    intent(in)   ::tr,r
DOUBLE PRECISION                  :: c_upper_fnc2(kdim), c_upper_temp(kdim)

INTEGER               :: i


do i=1,kdim

    c_upper_temp(i)=workcon(labor,capital,kprime(i),tr,r,effective_wage)
end do
c_upper_fnc2=c_upper_temp
END FUNCTION c_upper_fnc2

FUNCTION c_upper_fnc(prod,capital,kprime,kdim,tr,w,r)
IMPLICIT NONE
INTEGER, intent(in)   ::kdim
DOUBLE PRECISION,    intent(in)   ::kprime(kdim),capital
DOUBLE PRECISION,    intent(in)   ::tr,w,r,prod
DOUBLE PRECISION                  :: c_upper_fnc(kdim), c_upper_temp(kdim)

INTEGER               :: i

do i=1,kdim
    c_upper_temp(i)=workcon(0.90D0,capital,kprime(i),tr,r,w*prod)
end do
c_upper_fnc=c_upper_temp
END FUNCTION c_upper_fnc






FUNCTION workcon(labor,capital,kprime,tr,r,w)
IMPLICIT NONE
DOUBLE PRECISION,    intent(in)   :: capital,kprime,tr,r,labor,w
DOUBLE PRECISION                  :: workcon,ss_taxable,labor_tax,capital_tax

ss_taxable=min(ss_cap,labor*w)

workcon=(1+r)*(capital+tr)+labor*w-ss_taxable*tauss-kprime-all_taxer(r*(tr+capital),w*labor,1)+lump_sum

END FUNCTION workcon




FUNCTION rebate_calculator(labor_income,capital_income,young,labor_tax,capital_tax)
IMPLICIT NONE
INTEGER,    intent(in)   :: young
DOUBLE PRECISION,    intent(in)   :: labor_income,capital_income,labor_tax,capital_tax
DOUBLE PRECISION                  :: rebate_calculator,income,tax_comparison,ypp

ypp=labor_income-.5D0*tauss*min(labor_income,ss_cap)

if (young==0 .and. rebate_to_old == 0) then
    rebate_calculator=0.0D0
else
    if(labor_capital_rebate_base==1) then
        income=max(0.0D0,ypp)
    elseif(labor_capital_rebate_base==2) then
        income=max(0.0D0,capital_income)
    elseif(labor_capital_rebate_base==3) then
        income=max(0.0D0,ypp+capital_income)
    else
        print*,"DEFINE WHAT BASE FOR REBATE"
    END IF
    rebate_calculator=max(rebate_level+rebate_slope*income,0.0D0)
end if
END FUNCTION rebate_calculator






FUNCTION rebate_finder(labor_income,capital_income,young,labor_tax,capital_tax,rebate_level_test)
IMPLICIT NONE
INTEGER,    intent(in)   :: young
DOUBLE PRECISION,    intent(in)   :: labor_income,capital_income,labor_tax,capital_tax,rebate_level_test
DOUBLE PRECISION                  :: rebate_finder,income,tax_comparison,ypp

ypp=labor_income-.5D0*tauss*min(labor_income,ss_cap)

if (young==0 .and. rebate_to_old == 0) then
    rebate_finder=0.0D0
else
    if(labor_capital_rebate_base==1) then
        income=max(0.0D0,ypp)
        tax_comparison=labor_tax
    elseif(labor_capital_rebate_base==2) then
        income=max(0.0D0,capital_income)
        tax_comparison=capital_tax
    elseif(labor_capital_rebate_base==3) then
        income=max(0.0D0,ypp+capital_income)
        tax_comparison=labor_tax+capital_tax
    else
        print*,"DEFINE WHAT BASE FOR REBATE"
    END IF
    rebate_finder=max(rebate_level_test+rebate_slope*income,0.0D0)
end if

END FUNCTION rebate_finder



FUNCTION all_taxer(capital_income,labor_income,young)
IMPLICIT NONE
INTEGER,    intent(in)   :: young
DOUBLE PRECISION,    intent(in)   :: capital_income,labor_income
DOUBLE PRECISION                  :: all_taxer,ypp,rebate,capital_tax,labor_tax,old_labor_tax,tax_comparison
ypp=labor_income-.5D0*tauss*min(labor_income,ss_cap)

    labor_tax=lab_taxer(ypp)
    old_labor_tax=old_lab_taxer(ypp)
if(flat_tax==1) then
    if(ypp/avg_earnings_scale<lambda1) then
        labor_tax=0.0D0
    else
        labor_tax=old_lab_taxer(ypp)
    end if
end if
capital_tax=capital_income*lambdak
rebate=rebate_calculator(labor_income,capital_income,young,labor_tax,capital_tax)
if (no_increase==1 .and. labor_capital_rebate_base==1) then
    labor_tax=min(labor_tax-rebate,old_labor_tax)
else if(no_increase==-1) then
    labor_tax=labor_tax
else
    labor_tax=min(labor_tax,old_labor_tax)
end if

if(rebate_negative_tax==0) then
    if(labor_capital_rebate_base==1) then
        tax_comparison=labor_tax
    elseif(labor_capital_rebate_base==2) then
        tax_comparison=capital_tax
    elseif(labor_capital_rebate_base==3) then
        tax_comparison=labor_tax+capital_tax
    else
        print*,"DEFINE WHAT BASE FOR REBATE"
    END IF
    rebate=min(rebate,tax_comparison)
end if

all_taxer=labor_tax+capital_tax-rebate

end FUNCTION all_taxer


FUNCTION all_taxer_rebate(capital_income,labor_income,young,rebate_test)
IMPLICIT NONE
INTEGER,    intent(in)   :: young
DOUBLE PRECISION,    intent(in)   :: capital_income,labor_income,rebate_test
DOUBLE PRECISION                  :: all_taxer_rebate,ypp,rebate,capital_tax,labor_tax,old_labor_tax,tax_comparison
ypp=labor_income-.5D0*tauss*min(labor_income,ss_cap)
labor_tax=lab_taxer(ypp)
old_labor_tax=old_lab_taxer(ypp)
if(flat_tax==1) then
    if(ypp/avg_earnings_scale<lambda1) then
        labor_tax=0.0D0
    else
        labor_tax=old_lab_taxer(ypp)
    end if
end if
capital_tax=capital_income*lambdak
rebate=rebate_finder(labor_income,capital_income,young,labor_tax,capital_tax,rebate_test)

if (no_increase==1 .and. labor_capital_rebate_base==1) then
    labor_tax=min(labor_tax-rebate,old_labor_tax)
else if(no_increase==-1) then
    labor_tax=labor_tax
else
    labor_tax=min(labor_tax,old_labor_tax)
end if

if(rebate_negative_tax==0) then
    if(labor_capital_rebate_base==1) then
        tax_comparison=labor_tax
    elseif(labor_capital_rebate_base==2) then
        tax_comparison=capital_tax
    elseif(labor_capital_rebate_base==3) then
        tax_comparison=labor_tax+capital_tax
    else
        print*,"DEFINE WHAT BASE FOR REBATE"
    END IF
    rebate=min(rebate,tax_comparison)
end if

all_taxer_rebate=labor_tax+capital_tax-rebate

end FUNCTION all_taxer_rebate




FUNCTION all_taxer_b(capital_income,labor_income,young,lambda0_test,lambdak_test,rebate_test)
IMPLICIT NONE
INTEGER,    intent(in)   :: young
DOUBLE PRECISION,    intent(in)   :: capital_income,labor_income,lambda0_test,lambdak_test,rebate_test
DOUBLE PRECISION                  :: all_taxer_b,ypp,rebate,capital_tax,labor_tax,old_labor_tax,tax_comparison
ypp=labor_income-.5D0*tauss*min(labor_income,ss_cap)
    labor_tax=lab_taxer0(ypp,lambda0_test)
    old_labor_tax=old_lab_taxer(ypp)
if(flat_tax==1) then
    if(ypp/avg_earnings_scale<lambda1) then
        labor_tax=0.0D0
    else
        labor_tax=old_lab_taxer(ypp)
    end if
end if
capital_tax=capital_income*lambdak_test
rebate=rebate_finder(labor_income,capital_income,young,labor_tax,capital_tax,rebate_test)

if (no_increase==1 .and. labor_capital_rebate_base==1) then
    labor_tax=min(labor_tax-rebate,old_labor_tax)
else if(no_increase==-1) then
    labor_tax=labor_tax
else
    labor_tax=min(labor_tax,old_labor_tax)
end if

if(rebate_negative_tax==0) then
    if(labor_capital_rebate_base==1) then
        tax_comparison=labor_tax
    elseif(labor_capital_rebate_base==2) then
        tax_comparison=capital_tax
    elseif(labor_capital_rebate_base==3) then
        tax_comparison=labor_tax+capital_tax
    else
        print*,"DEFINE WHAT BASE FOR REBATE"
    END IF
    rebate=min(rebate,tax_comparison)
end if


all_taxer_b=labor_tax+capital_tax-rebate


end FUNCTION all_taxer_b

FUNCTION all_taxer_l(capital_income,labor_income,young,lambda0_test)
IMPLICIT NONE
INTEGER,    intent(in)   :: young
DOUBLE PRECISION,    intent(in)   :: capital_income,labor_income,lambda0_test
DOUBLE PRECISION                  :: all_taxer_l,ypp,rebate,capital_tax,labor_tax,old_labor_tax,tax_comparison
ypp=labor_income-.5D0*tauss*min(labor_income,ss_cap)
labor_tax=lab_taxer0(ypp,lambda0_test)
old_labor_tax=old_lab_taxer(ypp)
if(flat_tax==1) then
    if(ypp/avg_earnings_scale<lambda1) then
        labor_tax=0.0D0
    else
        labor_tax=old_lab_taxer(ypp)
    end if
end if

capital_tax=capital_income*lambdak
rebate=rebate_calculator(labor_income,capital_income,young,labor_tax,capital_tax)

if (no_increase==1 .and. labor_capital_rebate_base==1) then
    labor_tax=min(labor_tax-rebate,old_labor_tax)
else if(no_increase==-1) then
    labor_tax=labor_tax
else
    labor_tax=min(labor_tax,old_labor_tax)
end if

if(rebate_negative_tax==0) then
    if(labor_capital_rebate_base==1) then
        tax_comparison=labor_tax
    elseif(labor_capital_rebate_base==2) then
        tax_comparison=capital_tax
    elseif(labor_capital_rebate_base==3) then
        tax_comparison=labor_tax+capital_tax
    else
        print*,"DEFINE WHAT BASE FOR REBATE"
    END IF
    rebate=min(rebate,tax_comparison)
end if



all_taxer_l=labor_tax+capital_tax-rebate

end FUNCTION all_taxer_l

FUNCTION all_taxer_k(capital_income,labor_income,young,lambdak_test)
IMPLICIT NONE
INTEGER,    intent(in)   :: young
DOUBLE PRECISION,    intent(in)   :: capital_income,labor_income,lambdak_test
DOUBLE PRECISION                  :: all_taxer_k,ypp,rebate,capital_tax,labor_tax,old_labor_tax,tax_comparison
ypp=labor_income-.5D0*tauss*min(labor_income,ss_cap)
    labor_tax=lab_taxer(ypp)
    old_labor_tax=old_lab_taxer(ypp)
if(flat_tax==1) then
    if(ypp/avg_earnings_scale<lambda1) then
        labor_tax=0.0D0
    else
        labor_tax=old_lab_taxer(ypp)
    end if
end if
capital_tax=capital_income*lambdak_test
rebate=rebate_calculator(labor_income,capital_income,young,labor_tax,capital_tax)

if (no_increase==1 .and. labor_capital_rebate_base==1) then
    labor_tax=min(labor_tax-rebate,old_labor_tax)
else if(no_increase==-1) then
    labor_tax=labor_tax
else
    labor_tax=min(labor_tax,old_labor_tax)
end if

if(rebate_negative_tax==0) then
    if(labor_capital_rebate_base==1) then
        tax_comparison=labor_tax
    elseif(labor_capital_rebate_base==2) then
        tax_comparison=capital_tax
    elseif(labor_capital_rebate_base==3) then
        tax_comparison=labor_tax+capital_tax
    else
        print*,"DEFINE WHAT BASE FOR REBATE"
    END IF
    rebate=min(rebate,tax_comparison)
end if


all_taxer_k=labor_tax+capital_tax-rebate

end FUNCTION all_taxer_k


FUNCTION old_lab_taxer(income)
IMPLICIT NONE
DOUBLE PRECISION,    intent(in)   :: income
DOUBLE PRECISION                  :: old_lab_taxer,alt_taxer

if(income==0.0D0) then
    old_lab_taxer=0.0D0
else
    !lab_taxer=lambda0*(income-(income**(-lambda1)+lambda2)**(-1.0D0/lambda1))
old_lab_taxer=(1.0D0-lambda0_start*((income/avg_earnings_scale)**(-lambda1_start)))*income

end if
if(negative_taxes==0 .and. old_lab_taxer<0.0D0) then
    old_lab_taxer=0.0D0
end if

end FUNCTION old_lab_taxer


FUNCTION lab_taxer(income)
IMPLICIT NONE
DOUBLE PRECISION,    intent(in)   :: income
DOUBLE PRECISION                  :: lab_taxer,alt_taxer

if(income==0.0D0) then
    lab_taxer=0.0D0
else
    lab_taxer=(1.0D0-lambda0*((income/avg_earnings_scale)**(-lambda1)))*income
end if
if(negative_taxes==0 .and. lab_taxer<0.0D0) then
    lab_taxer=0.0D0
end if

end FUNCTION lab_taxer

FUNCTION lab_taxer0(income,lambda0_test)
IMPLICIT NONE
DOUBLE PRECISION,    intent(in)   :: income,lambda0_test
DOUBLE PRECISION                  :: lab_taxer0,alt_taxer

if(income==0.0D0) then
    lab_taxer0=0.0D0
else
    lab_taxer0=(1.0D0-lambda0_test*((income/avg_earnings_scale)**(-lambda1)))*income
end if
if(negative_taxes==0 .and. lab_taxer0<0.0D0) then
    lab_taxer0=0.0D0
end if

end FUNCTION lab_taxer0



FUNCTION consumption_tilda(con,eng)
IMPLICIT NONE
DOUBLE PRECISION,    intent(in)   :: con,eng
DOUBLE PRECISION                  :: consumption_tilda
consumption_tilda=con**gamma1*(eng-ebar)**(1.0D0-gamma1)
END FUNCTION consumption_tilda


FUNCTION retcon(capital,kprime,tr,r,ss_pass,avg_earnings)
IMPLICIT NONE
DOUBLE PRECISION,    intent(in)   :: capital,kprime,tr,r,ss_pass,avg_earnings
DOUBLE PRECISION                  :: retcon,tax_ss_income,labor_tax,capital_tax
    tax_ss_income=.0D0*ss_pass
    labor_tax=lab_taxer(tax_ss_income)
    capital_tax=lambdak*(r*(tr+capital))
    retcon=(1+r)*(capital+tr)+ss_pass-kprime+lump_sum-all_taxer(r*(tr+capital),tax_ss_income,0)
END FUNCTION retcon



FUNCTION retcon_ar(dim,capital,kprime,tr,r,ss_pass,avg_earnings)
IMPLICIT NONE
INTEGER,    intent(in)   :: dim
DOUBLE PRECISION,    intent(in)   :: capital,kprime(dim),tr,r,ss_pass,avg_earnings
DOUBLE PRECISION                  :: retcon_ar(dim),tax_ss_income,labor_tax,capital_tax
INTEGER               :: it
    tax_ss_income=.0D0*ss_pass

do it=1,dim
labor_tax=lab_taxer(tax_ss_income)
capital_tax=lambdak*(r*(tr+capital))

    retcon_ar(it)=(1+r)*(capital+tr)+ss_pass-kprime(it)+lump_sum-all_taxer(r*(tr+capital),tax_ss_income,0)
end do
END FUNCTION retcon_ar


FUNCTION utility_last(con,eng,ad_eqiv)
IMPLICIT NONE
DOUBLE PRECISION,    intent(in)   :: con,eng,ad_eqiv
DOUBLE PRECISION                  :: utils_tmp, utility_last,lc
if ( (eng+con)>0.015 .and. eng>ebar) then
    lc=con**gamma1*(eng-ebar)**(1.0D0-gamma1)/ad_eqiv
    utils_tmp =ad_eqiv*(lc**(1.0D0-sigma1)/(1.0D0-sigma1))
else
    utils_tmp = -1000000000000000000.0D0
end if
utility_last = utils_tmp
END FUNCTION utility_last

FUNCTION utility_last_cev(con,eng,cev,ad_eqiv,non_energy)
IMPLICIT NONE
INTEGER,    intent(in)   :: non_energy
DOUBLE PRECISION,    intent(in)   :: con,eng,ad_eqiv,cev
DOUBLE PRECISION                  :: utils_tmp, utility_last_cev,lc
if ( (eng+con)>0.015 .and. eng>ebar) then
    if(non_energy==1) then
        lc=con*((1.0D0+cev/100.0D0))**gamma1*(eng-ebar)**(1.0D0-gamma1)/ad_eqiv
    else
        lc=con**gamma1*(eng-ebar)**(1.0D0-gamma1)/ad_eqiv*(1.0D0+cev/100.0D0)
    end if
    utils_tmp =ad_eqiv*(lc**(1.0D0-sigma1)/(1.0D0-sigma1))
else
    utils_tmp = -1000000000000000000.0D0
end if
utility_last_cev = utils_tmp
END FUNCTION utility_last_cev


FUNCTION con_aggregator(con,eng)
IMPLICIT NONE
DOUBLE PRECISION,    intent(in)   :: con,eng
DOUBLE PRECISION                  :: lc,con_aggregator
if ( (eng+con)>0.015 .and. eng>ebar) then
    lc=con**gamma1*(eng-ebar)**(1.0D0-gamma1)
else
    lc = -1000000000000000000.0D0
end if
con_aggregator = lc
END FUNCTION con_aggregator


FUNCTION valuer_labor(labor,discount,vprime,ad_eqiv)
IMPLICIT NONE
DOUBLE PRECISION,    intent(in)   :: labor
DOUBLE PRECISION,    intent(in)   :: discount,vprime,ad_eqiv
DOUBLE PRECISION                  :: utils_tmp,valuer_labor
    if ( labor>0.01D0) then
        utils_tmp = 0.0D0-ad_eqiv*chi*labor**(1.0D0+1.0D0/sigma2)/(1.0D0+1.0D0/sigma2)+discount*beta*vprime
    else
        utils_tmp =0.0D0+discount*beta*vprime
    end if
    valuer_labor = utils_tmp
END FUNCTION valuer_labor

FUNCTION valuer(c,e,labor,discount,vprime,ad_eqiv)
IMPLICIT NONE
DOUBLE PRECISION,    intent(in)   :: c,e,labor,ad_eqiv
DOUBLE PRECISION,    intent(in)   :: discount,vprime
DOUBLE PRECISION                  :: utils_tmp,valuer,lc
    if ( (c+e)>0.015 .and. e>ebar) then
        lc=(c**gamma1*(e-ebar)**(1.0D0-gamma1))/ad_eqiv
        if ( labor>0.01D0) then
            utils_tmp = ad_eqiv*(lc**(1.0D0-sigma1)/(1.0D0-sigma1)-chi*labor**(1.0D0+1.0D0/sigma2)/(1.0D0+1.0D0/sigma2))+discount*beta*vprime
        else
            utils_tmp =ad_eqiv*(lc**(1.0D0-sigma1)/(1.0D0-sigma1))+discount*beta*vprime
        end if
    else
        utils_tmp = -1000000000000000000.0D0
    end if
    valuer = utils_tmp
END FUNCTION valuer

FUNCTION valuer_cev(c,e,labor,discount,vprime,ad_eqiv,cev,non_energy)
IMPLICIT NONE
INTEGER,    intent(in)   :: non_energy
DOUBLE PRECISION,    intent(in)   :: c,e,labor,cev,ad_eqiv
DOUBLE PRECISION,    intent(in)   :: discount,vprime
DOUBLE PRECISION                  :: utils_tmp,valuer_cev,lc
    if ( (c+e)>0.015 .and. e>ebar) then
        if(non_energy==1) then
            lc=((c*(1.0D0+cev/100.0D0))**gamma1*(e-ebar)**(1.0D0-gamma1))/ad_eqiv
        else
            lc=(c**gamma1*(e-ebar)**(1.0D0-gamma1))*(1.0D0+cev/100.0D0)/ad_eqiv
        end if
        if ( labor>0.01D0) then
            utils_tmp = ad_eqiv*(lc**(1.0D0-sigma1)/(1.0D0-sigma1)-chi*labor**(1.0D0+1.0D0/sigma2)/(1.0D0+1.0D0/sigma2))+discount*beta*vprime
        else
            utils_tmp =ad_eqiv*(lc**(1.0D0-sigma1)/(1.0D0-sigma1))+discount*beta*vprime
        end if
    else
        utils_tmp = -1000000000000000000.0D0
    end if
    valuer_cev = utils_tmp
END FUNCTION valuer_cev



FUNCTION valuer_ret(c,e,discount,vprime,ad_eqiv)
IMPLICIT NONE
DOUBLE PRECISION,    intent(in)   :: c,e,ad_eqiv
DOUBLE PRECISION,    intent(in)   :: discount
DOUBLE PRECISION                  :: utils_tmp, valuer_ret,vprime,lc
      if ( (c+e)>0.015 .and. e>ebar) then
        lc=c**gamma1*(e-ebar)**(1.0D0-gamma1)/ad_eqiv
        utils_tmp = ad_eqiv*lc**(1.0D0-sigma1)/(1.0D0-sigma1)+beta*discount*vprime
      else
         utils_tmp = -1000000000000000000.0D0
      end if
      valuer_ret = utils_tmp
END FUNCTION valuer_ret

FUNCTION valuer_ret_cev(c,e,discount,vprime,ad_eqiv,cev,non_energy)
IMPLICIT NONE
INTEGER,    intent(in)   :: non_energy
DOUBLE PRECISION,    intent(in)   :: c,e
DOUBLE PRECISION,    intent(in)   :: discount,cev,ad_eqiv
DOUBLE PRECISION                  :: utils_tmp, valuer_ret_cev,vprime,lc
    if ( (c+e)>0.015 .and. e>ebar) then
        if(non_energy==1) then
            lc=(c*(1.0D0+cev/100.0D0))**gamma1*(e-ebar)**(1.0D0-gamma1)/ad_eqiv
            utils_tmp = ad_eqiv*lc**(1.0D0-sigma1)/(1.0D0-sigma1)+beta*discount*vprime
        else
            lc=c**gamma1*(e-ebar)**(1.0D0-gamma1)*(1.0D0+cev/100.0D0)/ad_eqiv
            utils_tmp = ad_eqiv*lc**(1.0D0-sigma1)/(1.0D0-sigma1)+beta*discount*vprime
        end if
    else
        utils_tmp = -1000000000000000000.0D0
    end if
    valuer_ret_cev = utils_tmp
END FUNCTION valuer_ret_cev


SUBROUTINE energy_root_finder(inputs,misses,choicesN)
!
IMPLICIT NONE
INTEGER, intent(in)   ::choicesN
DOUBLE PRECISION,    intent(out)  :: misses(choicesN)
DOUBLE PRECISION,    intent(in)   :: inputs(choicesN)
DOUBLE PRECISION   :: Ep_guess,E_guess,K1_guess,K2_guess,N1_guess,N2_guess,K_in,N_in,Ec_in,Ae_in,Atfp_in,pe_in,taue_in
K_in=passed_parameters(1)
N_in=passed_parameters(2)
Ec_in=passed_parameters(3)
Atfp_in=passed_parameters(4)
Ae_in=passed_parameters(5)
taue_in=passed_parameters(6)

Ep_guess=inputs(1)
K1_guess=inputs(2)
N1_guess=inputs(3)

K2_guess=K_in-K1_guess
N2_guess=N_in-N1_guess
E_guess=Ae_in*K2_guess**alpha2*N2_guess**(1.0D0-alpha2)

pe_in=theta*Atfp_in*(K1_guess/Ep_guess)**alpha1*(N1_guess/Ep_guess)**(1.0D0-alpha1-theta)-taue_in
if(N1_guess<0.0D0 .or. N2_guess<0.0D0 .or. K1_guess<0.0D0 .or. K2_guess<0.0D0 .or. Ep_guess<0.0D0) then
misses(1)=100.0D0*(max(abs(N1_guess),0.0D0)+max(abs(N2_guess),0.0D0)+max(abs(K1_guess),0.0D0)+max(abs(K2_guess),0.0D0)+max(abs(E_guess),0.0D0)+max(abs(Ep_guess),0.0D0))
misses(2)=misses(1)*10.0D0
misses(3)=misses(1)*5.0D0
else
misses(1)=E_guess-Ep_guess-Ec_in
misses(2)=Atfp_in*alpha1*(N1_guess/K1_guess)**(1.0D0-alpha1-theta)*(Ep_guess/K1_guess)**theta-pe_in*alpha2*Ae_in*(N2_guess/K2_guess)**(1.0D0-alpha2)
misses(3)=Atfp_in*(1.0D0-alpha1-theta)*(K1_guess/N1_guess)**alpha1*(Ep_guess/N1_guess)**theta-pe_in*Ae_in*(1.0D0-alpha2)*(K2_guess/N2_guess)**alpha2
end if

END SUBROUTINE energy_root_finder



FUNCTION energy_share_equation(c_con,totconsumption,age)
IMPLICIT NONE
INTEGER, intent(in)   ::age
DOUBLE PRECISION,    intent(in)   :: c_con,totconsumption
DOUBLE PRECISION   :: e_con,energy_share_equation

e_con=(totconsumption-c_con)/(pe+taue)

energy_share_equation =&
-(((c_con)**gamma1*((e_con)-ebar)**(1.0D0-gamma1))**(1.0D0-sigma1)/(1.0D0-sigma1))

END FUNCTION energy_share_equation




DOUBLE PRECISION FUNCTION energy_sharer(ax,bx,cx,func,tol,xmin,totconsumption,age)
USE nrtype; USE nrutil, ONLY : nrerror
IMPLICIT NONE
INTEGER, intent(in)   ::age
    DOUBLE PRECISION, INTENT(IN)    :: ax,bx,cx,tol,totconsumption
    DOUBLE PRECISION, INTENT(OUT)   :: xmin
INTEGER, PARAMETER          :: ITMAX=100
DOUBLE PRECISION, PARAMETER ::CGOLD=0.3819660D0,ZEPS=1.0e-3_sp*epsilon(ax)
INTEGER                     ::iter
DOUBLE PRECISION            :: a,b,d,e,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm
DOUBLE PRECISION :: func
a=min(ax,cx)
b=max(ax,cx)
v=bx
w=v
x=v
e=0.0
fx=func(x,totconsumption,age)
fv=fx
fw=fx
do iter=1,ITMAX
    xm=0.5D0*(a+b)
    tol1=tol*abs(x)+ZEPS
    tol2=2.0D0*tol1
    if (abs(x-xm) <= (tol2-0.5D0*(b-a))) then
        xmin=x
        energy_sharer=fx
        RETURN
    end if
    if (abs(e) > tol1) then
        r=(x-w)*(fx-fv)
        q=(x-v)*(fx-fw)
        p=(x-v)*q-(x-w)*r
        q=2.0D0*(q-r)
        if (q > 0.0D0) p=-p
        q=abs(q)
        etemp=e
        e=d
        if (abs(p) >= abs(0.5D0*q*etemp) .or. &
            p <= q*(a-x) .or. p >= q*(b-x)) then
            e=merge(a-x,b-x, x >= xm )
            d=CGOLD*e
        else
            d=p/q
            u=x+d
            if (u-a < tol2 .or. b-u < tol2) d=sign(tol1,xm-x)
        end if
    else
        e=merge(a-x,b-x, x >= xm )
        d=CGOLD*e
    end if
    u=merge(x+d,x+sign(tol1,d), abs(d) >= tol1 )
    fu=func(u,totconsumption,age)
    if (fu <= fx) then
        if (u >= x) then
            a=x
        else
            b=x
        end if
        call shft(v,w,x,u)
        call shft(fv,fw,fx,fu)
    else
        if (u < x) then
            a=u
        else
            b=u
        end if
        if (fu <= fw .or. w == x) then
            v=w
            fv=fw
            w=u
            fw=fu
        else if (fu <= fv .or. v == x .or. v == w) then
            v=u
            fv=fu
        end if
    end if
end do
call nrerror('energy_sharer: exceed maximum iterations')
CONTAINS
!BL
SUBROUTINE shft(a,b,c,d)
    DOUBLE PRECISION, INTENT(OUT)   :: a
    DOUBLE PRECISION, INTENT(INOUT) :: b,c
    DOUBLE PRECISION, INTENT(IN)    :: d
a=b
b=c
c=d
END SUBROUTINE shft
END FUNCTION energy_sharer





DOUBLE PRECISION FUNCTION energy_leveler(ax,bx,cx,func,tol,xmin,K,N)
USE nrtype; USE nrutil, ONLY : nrerror
IMPLICIT NONE
    DOUBLE PRECISION, INTENT(IN)    :: ax,bx,cx,tol,K,N
    DOUBLE PRECISION, INTENT(OUT)   :: xmin
INTEGER, PARAMETER          :: ITMAX=100
DOUBLE PRECISION, PARAMETER ::CGOLD=0.3819660D0,ZEPS=1.0e-3_sp*epsilon(ax)
INTEGER                     ::iter
DOUBLE PRECISION            :: a,b,d,e,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm
DOUBLE PRECISION :: func
a=min(ax,cx)
b=max(ax,cx)
v=bx
w=v
x=v
e=0.0
fx=func(x,K,N)
fv=fx
fw=fx
do iter=1,ITMAX
    xm=0.5D0*(a+b)
    tol1=tol*abs(x)+ZEPS
    tol2=2.0D0*tol1
    if (abs(x-xm) <= (tol2-0.5D0*(b-a))) then
        xmin=x
        energy_leveler=fx
        RETURN
    end if
    if (abs(e) > tol1) then
        r=(x-w)*(fx-fv)
        q=(x-v)*(fx-fw)
        p=(x-v)*q-(x-w)*r
        q=2.0D0*(q-r)
        if (q > 0.0D0) p=-p
        q=abs(q)
        etemp=e
        e=d
        if (abs(p) >= abs(0.5D0*q*etemp) .or. &
            p <= q*(a-x) .or. p >= q*(b-x)) then
            e=merge(a-x,b-x, x >= xm )
            d=CGOLD*e
        else
            d=p/q
            u=x+d
            if (u-a < tol2 .or. b-u < tol2) d=sign(tol1,xm-x)
        end if
    else
        e=merge(a-x,b-x, x >= xm )
        d=CGOLD*e
    end if
    u=merge(x+d,x+sign(tol1,d), abs(d) >= tol1 )
    fu=func(u,K,N)
    if (fu <= fx) then
        if (u >= x) then
            a=x
        else
            b=x
        end if
        call shft(v,w,x,u)
        call shft(fv,fw,fx,fu)
    else
        if (u < x) then
            a=u
        else
            b=u
        end if
        if (fu <= fw .or. w == x) then
            v=w
            fv=fw
            w=u
            fw=fu
        else if (fu <= fv .or. v == x .or. v == w) then
            v=u
            fv=fu
        end if
    end if
end do
call nrerror('energy_leveler: exceed maximum iterations')
CONTAINS
!BL
SUBROUTINE shft(a,b,c,d)
    DOUBLE PRECISION, INTENT(OUT)   :: a
    DOUBLE PRECISION, INTENT(INOUT) :: b,c
    DOUBLE PRECISION, INTENT(IN)    :: d
a=b
b=c
c=d
END SUBROUTINE shft
END FUNCTION energy_leveler


subroutine spline_linear_val2( ndata, tdata, ydata, tval, yval, ndata2 )
      INTEGER :: left, right,i
      INTEGER, intent(in)  ::ndata,ndata2
      DOUBLE PRECISION,  intent(in) :: tdata(ndata),ydata(ndata),tval(ndata2)
      DOUBLE PRECISION,  intent(out) :: yval(ndata2)
      DOUBLE PRECISION  ::  ypval
ypval=0.0D0
yval=0.0D0
do i=1,ndata2
    call rvec_bracket2 ( ndata, tdata, tval(i), left, right )
  ypval = ( ydata(right) - ydata(left) ) / ( tdata(right) - tdata(left) )
  yval(i) = ydata(left) +  ( tval(i) - tdata(left) ) * ypval
end do
  return
end SUBROUTINE spline_linear_val2

subroutine rvec_bracket2 ( n, x, xval, left, right )
  implicit none
      INTEGER, intent(in)  ::n
      INTEGER, intent(out)  ::left,right
      DOUBLE PRECISION, intent(in)   ::xval,x(n)
  integer i
  do i = 2, n - 1

    if ( xval < x(i) ) then
      left = i - 1
      right = i
      return
    end if
   end do
  left = n - 1
  right = n
  return
end subroutine rvec_bracket2

END module functions_mod

